#include <sys/time.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "compute.h"


void do_compute(const struct parameters* p, struct results *r)
{
 
	// put them on the stack so it will be faster.
	const size_t N = p->N;
	const size_t M = p->M; 

	//temperature matrix
	double (*t)[N+2][M+2];
	//conductivity matrix
	double (*c)[N][M];

	t = (double (*)[N+2][M+2]) p->tinit;
	c = (double (*)[N][M]) p->conductivity;
	
	size_t i;
	size_t j;
	
	//row used to save the values calculated until they can be written. 
	double temp_row[M+1];
	// auxiliar variables used to save temporary results
	double aux, auxdiff;
	
	// weights of the direct neighbors.
	const double direct_neigh_weight = 0.146446609;  // sqrt(2)/(sqrt(2) +1) /4;
	const double diag_neigh_weight = 0.103553391; // 1/(sqrt(2) +1) /4;
	//starting average value of the values from the matrix
	double avg = (p->io_tmax + p->io_tmin) * 0.5;
	//sum and average sum of differences between avg and the values in 1 row of the matrix.  
	double delta, avg_delta;	
	
	const double inverseM = 1.0/M;
	const double inverseN = 1.0/N;
	//max difference between old and new values of the temperatures.
	double maxdiff;

	struct timeval tim;
	//start time at the beginning of the loop
    double loop_start; 
	//end time at the ending of the loop
    double loop_end;
	double prog_start, prog_end;
	
	gettimeofday(&tim, NULL);
	prog_start=tim.tv_sec+(tim.tv_usec/1000000.0);


/******************************************************************************************************************/
/*
 for (i = 0; i< N; i++)
 {
 	for(j = 0; j < M; j++)
 	{
 		printf("%lf ", (*c)[i][j]);
 	}
 	printf("\n");
 }
 */
 
/******************************************************************************************************************/
 for (i = 0; i<= N+1; i++)
	{
	for (j = 0 ; j<= M+1; j++)
	{
	printf( "%lf  ",(*t)[i][j]);
	if( j == 0 || j == M)
	printf (" * "); 
 	}
 	printf("\n");
	if( i == 0 || i ==  N  )
	 printf("\n");
 } 
 
	
	do	{
		gettimeofday(&tim, NULL);
		loop_start=tim.tv_sec+(tim.tv_usec/1000000.0);
		//re-init variables
		r->tmin = p->io_tmax;
		r->tmax = p->io_tmin;
		delta = 0;
		maxdiff = 0;

		//calculate the new values for the first row. 
		for( j = 1; j <= M; ++j)
		{
			 temp_row[j] = (*c)[0][j-1] *  (*t)[1][j] + 
							(1 - (*c)[0][j-1] ) * 
				 					( 
			 						direct_neigh_weight * ( (*t)[1][j-1] + (*t)[1][j+1] + (*t)[0][j] + (*t)[2][j] ) +
				 					diag_neigh_weight * ( (*t)[0][j+1] + (*t)[0][j-1] + (*t)[2][j+1] + (*t)[2][j-1] )
				 					);

		 	 delta = delta + temp_row[j] - avg; 
		 	 
		 	 auxdiff = fabs((*t)[1][j] -  temp_row[j] );
		 	 if(auxdiff > maxdiff)
		 	 	maxdiff = auxdiff;
			 
			 if(r->tmin > temp_row[j] ) 
			 {
				r->tmin = temp_row[j];
	
			 }else if(r->tmax < temp_row[j] )
			 {
			  	r->tmax = temp_row[j];
		 	 }
		}
		
		//important step. Initialize the first parameter with the left diagonal value of the next point  
		temp_row[0] = (*t)[1][0];
		//calculate average delta for the row. 
		avg_delta = delta * inverseM;
		//at the first iteration update the avg to the avg of the first row.
		if(r->niter ==0){
			avg = avg + avg_delta;
			avg_delta = 0;
		}
		//calculate the new values for the rest of the matrix. 
		for( i = 2; i <= N; ++i)
		{
			delta = 0;	
			for( j = 1; j <= M; ++j)
			{
				aux = 	(*c)[i-1][j-1] * (*t)[i][j] + 
						(1 -(*c)[i-1][j-1] ) * 
									( 
							 		direct_neigh_weight * ( (*t)[i][j-1] + (*t)[i][j+1] + (*t)[i-1][j] + (*t)[i+1][j] ) +
							 		diag_neigh_weight * ( (*t)[i-1][j+1] + temp_row[0] + (*t)[i+1][j+1] + (*t)[i+1][j-1] )
							 		);
				temp_row[0] = (*t)[i-1][j];
				(*t)[i-1][j] = temp_row[j];
				temp_row[j] = aux;
				delta = delta + aux - avg;

				auxdiff = fabs((*t)[i][j] -  aux);
		 	 	if(auxdiff > maxdiff)
			 	 	maxdiff = auxdiff;
				
				if(r->tmin > aux ) 
				{
					r->tmin = aux;
	
				}else if(r->tmax < aux )
				{
					r->tmax = aux;
				}
			}
			avg_delta = avg_delta + delta * inverseM;
			// prepare the diagonal value for the first element from the next row
			temp_row[0] = (*t)[i][0]; 
		}
		//copy the values calculated for the last row.
		for(j = 1; j <= M; ++j)
			(*t)[N][j]  = temp_row[j];
		
		//update the borders of the matrix ( the values on the border columns).	
		for( i =1 ; i <=N; ++i)	
		{
			(*t)[i][0] = (*t)[i][M];
			(*t)[i][M+1] = (*t)[i][1];
		}

		//update reported values
		r->niter = r->niter + 1;
		avg = avg + inverseN * avg_delta; 
		r->tavg = avg;
		r->maxdiff = maxdiff;
		gettimeofday(&tim, NULL);
   		loop_end=tim.tv_sec+(tim.tv_usec/1000000.0);
 		r->time =  loop_end - loop_start;
		
		report_results(p,r);
		}while(r->niter < p->maxiter && p->threshold < r->maxdiff );
//end of algorithm.
	
	gettimeofday(&tim, NULL);
	prog_end =tim.tv_sec+(tim.tv_usec/1000000.0);
	
	printf("\n TOTAL program execution time  : %lf\n", prog_end - prog_start);
/******************************************************************************************************************/
/*	printf("\n");
	for (i = 0; i<= N+1; i++)
	{
		for (j = 0 ; j<= M+1; j++)
		{
			printf( "%lf  ",(*t)[i][j]);
			if( j == 0 || j == M)
				printf (" * "); 
	 	}
	 		printf("\n");
		if( i == 0 || i ==  N  )
			 printf("\n");
	 } 
*/ }

